# ----------------------------------------------------------------------
# How Does Shaming Human Rights Violators Abroad Shape Attitudes at Home?
# Study I: Script to produce figures for manuscript and SM
# ----------------------------------------------------------------------

# ----------------------------------------------------------------------
# load all relevant themes and packages
# ----------------------------------------------------------------------
rm(list=ls())
library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("RItools")
library("estimatr")
library("fastDummies")
library("ggpubr")
library("estimatr")
library("mediation")
library("car")
library("table1")
library("psy")
library("dplyr")
library("ggstatsplot")
library("ggplot2")
library("stringr")
library("tidyverse")
library("haven")
library("lfe")
library("reshape2")
library("naniar")
library("fixest")
library("Hmisc")
library("quanteda")
library("readtext")
library("quanteda.textplots")
library("htmltools")
library("devtools")
install_github("naoki-egami/exr", dependencies = TRUE)
library("exr")

# ----------------------------------------------------------------------
# Read data
# ----------------------------------------------------------------------

data <- read_csv("Data_S1.csv")
data <- data[3:nrow(data),]

# ----------------------------------------------------------------------
# Re-code covariates and outcomes
# ----------------------------------------------------------------------

#OUTCONES summary and rescaling

#`gov_approve` -- rescale: high # indicates support for government
data$gov_approve = recode(data$gov_approve1, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

#`torture` -- high # indicate less support for torture (no need to rescale)

# `legal` -- index: create index of legal obligation measure (Bayram 2017)
data <- data %>% 
  mutate(.,
         legal = (as.numeric(data$legal_obligation_1)+
                    as.numeric(data$legal_obligation_2)+
                    as.numeric(data$legal_obligation_3)+
                    as.numeric(data$legal_obligation_4)+
                    as.numeric(data$legal_obligation_5))/5)

#rescale so that higher #s mean more support
data$legal<-6-data$legal

# `HR values` -- high # indicates more respect for human rights (no change required)

# morality -- rescale: high # indicates U.S. more moral 
data$morality = recode(data$morality, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

# power -- rescale: high # indicates U.S. more powerful 
data$power = recode(data$power, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

#create a Hawk and Dove index with items 1,2,4,5 (3 was an attention check)
#rescale item 2 to be in the same direction as others (more hawkish)
data$attention2 = recode(data$attention2, "1=5; 2=4; 3=3; 4=2; 5=1")
data <- data %>% 
  mutate(.,hawk= (as.numeric(data$attention1)+
            as.numeric(data$attention2)+
            as.numeric(data$attention4)+
            as.numeric(data$attention5))/4)
data$partisanship
        
# ----------------------------------------------------------------------
# recode variables
# ----------------------------------------------------------------------
data <- data %>%
  mutate(.,
         treatment = case_when(
           condition == "no_shame" ~ 0,
           condition == "shame" ~ 1),
         Democrat = case_when(
           partisanship == 1 ~ 1,
           partisanship == 2 ~ 0,
           other == 3 ~ 0,
           other == 1 ~ 0,
           other == 2 ~ 1),
         Republican = case_when(
           partisanship == 1 ~ 0,
           partisanship == 2 ~ 1,
           other == 3 ~ 0,
           other == 1 ~ 1,
           other == 2 ~ 0),
         Independents = case_when(
           other == 3 ~ 1,
           partisanship == 1 ~ 0,
           partisanship == 2 ~ 0,
           other == 1 ~ 0,
           other == 2 ~ 0
         ),
         Female = ifelse(gender == 2, 1, 0),
         White = ifelse(race == 1, 1, 0),
         Black = ifelse(race == 2,1,0),
         Hispanic = ifelse(race==3,1,0),
         Asian = ifelse(race==4,1,0),
         Native_american = ifelse(race==5,1,0),
         Middel_eastern = ifelse(race==6,1,0),
         Mixed = ifelse(race==7,1,0))

# ----------------------------------------------------------------------
# Produce Figure 4 (manuscript)
# ----------------------------------------------------------------------

#H1: does shaming affect government approval?
H1<- lm_robust(gov_approve~treatment, data = data)%>% tidy %>% 
  mutate(.,
         model = "Government \n Approval \n (H1)",
         Outcome = "a") 


#H2: does shaming affect human rights perception?
H2<- lm_robust(`HR values`~treatment, data=data)%>% tidy %>% 
  mutate(.,
         model = "Human Rights \n Perception \n (H2)",
         Outcome = "a")


#H3a: does shaming affect opposition of torture?
H3<- lm_robust(torture~treatment, data = data) %>% tidy %>% 
  mutate(.,
         model = "Oppose \n Torture \n (H3)",
         Outcome = "a")


#Create figure for main effects

main_fx <- rbind(H1, H2, H3) %>% 
  filter(.,
         term == "treatment")

order_x <- c("Government \n Approval \n (H1)", "Human Rights \n Perception \n (H2)",
             "Oppose \n Torture \n (H3)" 
             )
# Generate coefficient plot 
ggplot(main_fx, aes(x = model, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4), color = 'darkblue',
                  fill = "white") +
  geom_text(aes(label= round(estimate, digits = 1)), hjust=-.6, size = 3, color = "firebrick3",
            family = "Times") +
  geom_text(aes(label= paste("(",round(std.error, digits = 1),")")), hjust=-.2, vjust =2 ,  size = 3, 
            family = "Times") +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("../Main_Figs/fig4.pdf", width = 6, height = 4)

#substantive discussion of Fig 4:
#what was average gov approval in each condition?
data%>%filter(treatment=="0")%>%{mean(as.numeric(.$gov_approve),na.rm=T)} #3
data%>%filter(treatment=="1")%>%{mean(as.numeric(.$gov_approve),na.rm=T)} #5.3
#effect in terms of SDs?
2.3/sd(data$gov_approve,na.rm = T)

#what was average HR perceptions in each condition?
data%>%filter(treatment=="0")%>%{mean(as.numeric(.$`HR values`),na.rm=T)} #3.7
data%>%filter(treatment=="1")%>%{mean(as.numeric(.$`HR values`),na.rm=T)} #5.4
#effect in terms of SDs?
1.7/sd(data$`HR values`,na.rm = T)

#what was average torture in each condition?
data%>%filter(treatment=="0")%>%{mean(as.numeric(.$torture),na.rm=T)} #5.8
data%>%filter(treatment=="1")%>%{mean(as.numeric(.$torture),na.rm=T)} #5.6
#effect in terms of SDs?
0.2/sd(data$torture,na.rm = T)

# ----------------------------------------------------------------------
# Produce Figure 5 (manuscript)
# ----------------------------------------------------------------------

#Perception of Morality
morality<- lm_robust(morality~treatment, data = data) %>% tidy %>% 
  mutate(.,
         model = "Perception of \n U.S. Morality",
         Outcome = "a")

#Perception of Power
power<- lm_robust(power~treatment, data = data) %>% tidy %>% 
  mutate(.,
         model = "Perception of \n U.S. Power",
         Outcome = "a")

#Create figure

main_fx <- rbind(morality, power) %>% 
  filter(.,
         term == "treatment")

order_x <- c("Perception of \n U.S. Morality", "Perception of \n U.S. Power")
# Generate coefficient plot 
ggplot(main_fx, aes(x = model, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4), color = 'darkblue',
                  fill = "white") +
  geom_text(aes(label= round(estimate, digits = 1)), hjust=-.6, size = 3, color = "firebrick3",
            family = "Times") +
  geom_text(aes(label= paste("(",round(std.error, digits = 1),")")), hjust=-.2, vjust =2 ,  size = 3, 
            family = "Times") +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("../Main_Figs/fig5.pdf", width = 6, height = 4)

#substantive discussion of Fig 5:
#what was average morality peception in each condition?
data%>%filter(treatment=="0")%>%{mean(as.numeric(.$morality),na.rm=T)} #3.5
data%>%filter(treatment=="1")%>%{mean(as.numeric(.$morality),na.rm=T)} #5.3
#effect in terms of SDs?
1.8/sd(data$morality,na.rm = T)

#what was average power peception in each condition?
data%>%filter(treatment=="0")%>%{mean(as.numeric(.$power),na.rm=T)} #3.8
data%>%filter(treatment=="1")%>%{mean(as.numeric(.$power),na.rm=T)} #4.8
#effect in terms of SDs?
1.1/sd(data$power,na.rm = T)

# ----------------------------------------------------------------------
# Supplementary Materials
# ----------------------------------------------------------------------

#All files are saved into the Appendix/S1 folders

# ----------------------------------------------------------------------
# Descriptive statistics
# ----------------------------------------------------------------------

data1 <- data %>% 
  filter(.,
         !is.na(partisanship))

data1$age2<-as.numeric(data1$age2)
data1$gov_approve<-as.numeric(data1$gov_approve)
data1$`HR values`<-as.numeric(data1$`HR values`)
data1$torture<-as.numeric(data1$torture)
data1$legal<-as.numeric(data1$legal)
data1$power<-as.numeric(data1$power)
data1$morality<-as.numeric(data1$morality)

disc_table <- dplyr::select(data1, c(gov_approve,`HR values`, torture, legal, morality, power,
                                    Democrat, Independents, Republican, Female, White, Black,
                                    Hispanic, Asian, Native_american, Middel_eastern, Mixed, age2))

stargazer(as.data.frame(disc_table),
          covariate.labels = c("Government Approval","Human Rights Respect","Oppose Torture", "Legal Obligation", "Morality", "Power",
                                "Democrat", "Independent", "Republican", "Female", 
                               "White", "Black/African American", "Hispanic", "Asian/Asian American", "Native American",
                               "Middle Eastern", "Mixed Race","Age"),
          
          title = "Descriptive Statistics - Survey Experiment I",
          label = "tab:dsc_exp",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="../Appendix/Tables/S1/A1.tex")

# ----------------------------------------------------------------------
# Table for main models
# ----------------------------------------------------------------------

#Create table of main models and controls
#Add (pre-registered) legal obligation outcome

H1<- lm(gov_approve~treatment, data = data)#approval
H2 <- lm(`HR values`~treatment, data = data)#perceptions
H3a<- lm(torture~treatment, data = data)#torture
H3b<- lm(legal~treatment, data = data)#obligation

##add demographic controls

#Govermnemt approval
H1_control<- lm(gov_approve~treatment + as.factor(gender) + as.factor(age2) 
                + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                  as.factor(political_party) + as.factor(region), data = data) 
summary(H1_control)

#HR perception
H2_control <- lm(`HR values`~ treatment + as.factor(gender) + as.factor(age2) 
                 + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                   as.factor(political_party) + as.factor(region), data = data) 
summary(H2_control)

#Torture
H3a_control<- lm(torture~treatment + as.factor(gender) + as.factor(age2) 
                 + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                   as.factor(political_party) + as.factor(region), data = data) 
summary(H3a_control)

#Legal
H3b_control<- lm(legal~treatment + as.factor(gender) + as.factor(age2) 
                 + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                   as.factor(political_party) + as.factor(region), data = data) 
summary(H3b_control)


##Create table with controls

#Create table
tab1<-stargazer(H1, H1_control, H2, H2_control, H3a, H3a_control, H3b, H3b_control,
                out= "../Appendix/Tables/S1/A2.tex",  style="qje",
                title="Main Models with and without (pre-treatment) Demographic Controls",
                keep = "treatment",
                dep.var.labels = c("Gov Approval", "HR Respect", "Oppose Torture", "Legal Obligation"),
                covariate.labels = "Treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:controls",
                # star.char = c("", "", ""),
                add.lines = list(
                  c("Controls", "No", "Yes","No", "Yes", "No","Yes", "No", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Demographic controls"))

note.latex <- "\\multicolumn{9}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S1/A2.tex')


# ----------------------------------------------------------------------
# Mediation analysis - exploratory
# ----------------------------------------------------------------------

##Fig A1 (mediator: morality, outcome: torture)##
H3_med <- lm(morality~treatment, 
             data=data)
summary(H3_med)
data$torture<-as.numeric(data$torture)
H3_out <- lm(torture~treatment + morality, data = data)
summary(H3_out)

mediation_t <- mediate(H3_med, H3_out, treat = "treatment", 
                              mediator = "morality", dropobs=TRUE, 
                              boot=TRUE, conf.level=.90)

summary(mediation_t)

#save pdf
pdf(file="../Appendix/Figures/S1/A1.pdf")
plot(mediation_t, labels=c("ACME\n(Morality)",
                            "Direct \nEffect \n(Torture)", "Total \nEffect"), xlim = c(-2,1))
dev.off()


##Sensitivity analysis (Figure A5)
sens.out <- medsens(mediation_t, rho.by = 0.1, effect.type = "indirect", sims = 100) 

par(mar=c(5,6.5,3.2,01))

pdf(file="../Appendix/Figures/S1/A5.pdf", width = 13)
par( mfrow= c(1,2))
plot(sens.out, sign.prod="positive", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(a)")

plot(sens.out, sign.prod="negative", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(b)")
dev.off()


##Fig A2a (mediator: morality, outcome: approval)##

morality_med <- lm(morality~treatment, 
                     data=data)
summary(morality_med)

morality_out <- lm(gov_approve~treatment + morality, data = data)

mediation_morality <- mediate(morality_med, morality_out, treat = "treatment", 
                                mediator = "morality", dropobs=TRUE, 
                                boot=TRUE, conf.level=.90)
summary(mediation_morality)


#save pdf
pdf(file="../Appendix/Figures/S1/A2a.pdf")
plot(mediation_morality, labels=c("ACME\n(Morality)",
                                  "Direct \nEffect", "Total \nEffect"), xlim = c(0,3))
dev.off()

##Sensitivity analysis (A3)
sens.out <- medsens(mediation_morality, rho.by = 0.1, effect.type = "indirect", sims = 100) 
par(mar=c(5,6.5,3.2,01))

pdf(file="../Appendix/Figures/S1/A3.pdf", width = 15)
par( mfrow= c(1,2))
plot(sens.out, sign.prod="positive", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(a)")

plot(sens.out, sign.prod="negative", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(b)")
dev.off()


##Fig A2b (mediator: power, outcome: approval)##

##Mediation analysis
power_med <- lm(power~treatment, 
                   data=data)
summary(power_med)

power_out <- lm(gov_approve~treatment + power , data = data)

mediation_power <- mediate(power_med, power_out, treat = "treatment", 
                              mediator = "power", dropobs=TRUE, 
                              boot=TRUE, conf.level=.90)
summary(mediation_power)

#save pdf
pdf(file="../Appendix/Figures/S1/A2b.pdf")
plot(mediation_power, labels=c("ACME\n(Power)",
                                  "Direct \nEffect", "Total \nEffect"), xlim = c(0,3))
dev.off()

##Sensitivity analysis (A4)
sens.out <- medsens(mediation_power, rho.by = 0.1, effect.type = "indirect", sims = 100) 
par(mar=c(5,6.5,3.2,01))

pdf(file="../Appendix/Figures/S1/A4.pdf", width = 15)
par( mfrow= c(1,2))
plot(sens.out, sign.prod="positive", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(a)")

plot(sens.out, sign.prod="negative", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(b)")
dev.off()


# ----------------------------------------------------------------------
# Manipulation Checks
# ----------------------------------------------------------------------

#Did the U.S. government criticize the other country?

data$manipulation = recode(data$attention_2, "1=3; 2=2; 3=1") #recode such that [yes(3), not specified (2), no (1)]

#recode such that [yes(3), not specified (2), no (1)]

manipulation<- lm(manipulation~treatment, data = data) 
summary(manipulation)

#Respondents who got treatment are more likely to say the U.S. government criticized the other country.
manipulation_controls<-lm(manipulation~treatment+ as.factor(gender) + as.factor(age2) 
                          + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                            as.factor(political_party) + as.factor(region),data=data)

summary(manipulation_controls) #similar result when controlling for demographics

tab1<-stargazer(manipulation, manipulation_controls,
                out= "../Appendix/Tables/S1/A3.tex",  style="qje",
                title="Manipulation Check",
                keep = "treatment",
                dep.var.labels = "Belief that the U.S. government criticized the other country",
                covariate.labels = "Treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:manipulation",
               # star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
               omit.table.layout = NULL,
               notes.append = FALSE, 
               notes.align = "l",
               notes = c("Demographic controls"))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S1/A3.tex')

#how many answered manipulation check correctly?
data<-data%>%
  mutate(.,
         pass = case_when(attention_2=="1" & treatment == "1"~"1",
                          attention_2=="2" & treatment == "1"~"0",
                          attention_2=="3" & treatment == "1"~"0",
                          attention_2=="1" & treatment == "0"~"0",
                          attention_2=="2" & treatment == "0"~"0",
                          attention_2=="3" & treatment == "0"~"1"))

sum(as.numeric(data$pass),na.rm=T)/sum(!is.na(data$pass))

# ----------------------------------------------------------------------
# Exploring confounding
# ----------------------------------------------------------------------
#how many people thought of a specific country?
data1$placebo = ifelse(!is.na(data1$placebo_1_TEXT),1,0)
sum(data1$placebo)/nrow(data1) #617 -- about 45%

#Were people in the treatment group more likely to think of a specific country?
counfoudning<-lm(placebo~treatment, data=data1)
summary(counfoudning) #no

counfoudning_controls<-lm(placebo~treatment+ as.factor(gender) + as.factor(age2) 
                     + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                       as.factor(political_party) + as.factor(region),data=data1)

summary(counfoudning_controls) #similar result when controling for demographics

#Create table
tab1<-stargazer(counfoudning, counfoudning_controls,
                out= "../Appendix/Tables/S1/A4.tex",  style="qje",
                title="The effect of shaming treatment on thinking of a specific country",
                keep = "treatment",
                dep.var.labels = "Did you think of a specific country?",
                covariate.labels = "Treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:confounding",
                #star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Demographic controls"))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S1/A4.tex')

#Who did they think about? did they all think of the same country (Russia)?
placebo<-data[which(!is.na(data$placebo_1_TEXT)), ]
myvars <- "placebo_1_TEXT"
placebo <- placebo[myvars]
placebo$placebo_1_TEXT<-as.character(placebo$placebo_1_TEXT)
a<- corpus(placebo$placebo_1_TEXT)
b <- dfm(a, remove = stopwords("english"), remove_punct = TRUE)

topfeatures(b, 10) # 10 most frequent words

#Save wordcloud
png(file="../Appendix/Figures/S1/A6.png", width = 500, height = 500)
set.seed(100)
textplot_wordcloud(b, min_count = 10, random_order = FALSE,
                   rotation = .5, min_size = 2,
                   max_size = 8,
                   color = RColorBrewer::brewer.pal(8, "Dark2"))
dev.off()


#plot country response by the outcome question (Figure A7) 
#treatment effect by Russia, China, all other countries, no country

library(tidytext)
russia_tx <- data %>%
  filter(placebo_1_TEXT=="russia" | placebo_1_TEXT=="Russia" |
           placebo_1_TEXT=="RUSSIA")

china_tx <- data %>%
  filter(placebo_1_TEXT=="China" | placebo_1_TEXT=="china"|
           placebo_1_TEXT=="CHINA")  

other_tx <- data %>%
  filter(placebo_1_TEXT!="China" & placebo_1_TEXT!="china" & placebo_1_TEXT!="CHINA" &
           placebo_1_TEXT!="russia" & placebo_1_TEXT!="Russia" & placebo_1_TEXT!="RUSSIA") 

no_tx <- data %>%
  filter(is.na(placebo_1_TEXT))

#plot treatment effect by these 4 categories

#H1s: gov_approve~treatment
H1_russia<- lm_robust(gov_approve~treatment, data = russia_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Russia \n N=183",
         Outcome = "Gov Approval (H1)") 

H1_china<- lm_robust(gov_approve~treatment, data = china_tx)%>% tidy %>% 
  mutate(.,
         Sample = "China \n N=174",
         Outcome = "Gov Approval (H1)") 

H1_other<- lm_robust(gov_approve~treatment, data = other_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Other \n N=260",
         Outcome = "Gov Approval (H1)") 

H1_none<- lm_robust(gov_approve~treatment, data = no_tx)%>% tidy %>% 
  mutate(.,
         Sample = "None \n N=735",
         Outcome = "Gov Approval (H1)") 

#H2: `HR values`~treatment
H2_russia<- lm_robust(`HR values`~treatment, data=russia_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Russia \n N=183",
         Outcome = "HR Perception (H2)")

H2_china<- lm_robust(`HR values`~treatment, data=china_tx)%>% tidy %>% 
  mutate(.,
         Sample = "China \n N=174",
         Outcome = "HR Perception (H2)")

H2_other<- lm_robust(`HR values`~treatment, data=other_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Other \n N=260",
         Outcome = "HR Perception (H2)")

H2_none<- lm_robust(`HR values`~treatment, data=no_tx)%>% tidy %>% 
  mutate(.,
         Sample = "None \n N=735",
         Outcome = "HR Perception (H2)")

#H3a: torture~treatment
H3_russia<- lm_robust(torture~treatment, data=russia_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Russia \n N=183",
         Outcome = "Oppose torture (H3)")

H3_china<- lm_robust(torture~treatment, data=china_tx)%>% tidy %>% 
  mutate(.,
         Sample = "China \n N=174",
         Outcome = "Oppose torture (H3)")

H3_other<- lm_robust(torture~treatment, data=other_tx)%>% tidy %>% 
  mutate(.,
         Sample = "Other \n N=260",
         Outcome = "Oppose torture (H3)")

H3_none<- lm_robust(torture~treatment, data=no_tx)%>% tidy %>% 
  mutate(.,
         Sample = "None \n N=735",
         Outcome = "Oppose torture (H3)")

#Create figure for main effects
main_fx <- rbind(H1_russia, H1_china, H1_other,H1_none,
                 H2_russia, H2_china, H2_other,H2_none,
                 H3_russia, H3_china, H3_other,H3_none
                 ) %>% 
  filter(.,
         term == "treatment")

order_x <- c("Russia \n N=183", "China \n N=174", "Other \n N=260", "None \n N=735")

# Generate coefficient plot 
ggplot(main_fx, aes(x = Sample, y = estimate, color = Outcome, shape = Outcome)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4)) +
  labs(x = "Country Reported in Open-ended Question",
       y = "Treatment Effect Size") +
  scale_x_discrete(limits= order_x) +
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("../Appendix/Figures/S1/A7.pdf", width = 6, height = 4)


# ----------------------------------------------------------------------
# Attrition
# ----------------------------------------------------------------------

#how many attrited?
#data1 indicates all respondents who passed attention. Out of these -- how many attrited?
data1 <- data1 %>%
  mutate(., attrite = ifelse(is.na(data1$religion),1,0)) #religion is the last question in the survey

sum(data1$attrite) #62 respondents
(sum(data1$attrite)/nrow(data1))*100 #only 4.5% of the sample

#Does my treatment cause people to attrite?
attrite_model<-lm(attrite~treatment,data=data1)
summary(attrite_model) #no

attrite_controls<-lm(attrite~treatment+ as.factor(gender) + as.factor(age2) 
                + as.factor(ethnicity) + as.factor(education2) + as.factor(hispanic) +
                  as.factor(political_party) + as.factor(region),data=data1)
summary(attrite_controls) #also when controlling for demographics

#Create table
tab1<-stargazer(attrite_model, attrite_controls,
                out= "../Appendix/Tables/S1/A5.tex",  style="qje",
                title="The effect of shaming treatment on attrition",
                keep = "treatment",
                dep.var.labels = "Attrition",
                covariate.labels = "Treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:attrite",
               # star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Demographic controls"))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S1/A5.tex')

# ----------------------------------------------------------------------
# Probing External Validity
# ----------------------------------------------------------------------

#Create figure C11
#implement external validity bias using Devaux and Egami

#rename
data$HR<-data$`HR values`

covariates <- c("Democrat","Republican","Independents","sex", "White", "Black",
                "Hispanic", "Asian", "Native_american", "Middel_eastern", "Mixed", "age",
                "education","hawk"
) 

for_sate <- as.formula(paste0("HR ~ treatment + ", paste(covariates, collapse = "+")))
lm_sate <- lm(for_sate, data = data)
sate_est <- summary(lm_sate)$coef["treatment", c(1,2)]
sate_est

exr_out <- exr(outcome = "HR", 
               treatment = "treatment", 
               covariates = c("Democrat","Republican","Independents", "sex", "White", "Black",
                              "Hispanic", "Asian", "Native_american", "Middel_eastern", "Mixed", "age",
                              "education","hawk"), 
               data = data,
               boot = TRUE,
               sate_estimate = sate_est) 


plot(exr_out)

pdf(file="../Appendix/Figures/S1/C11b.pdf")
set.seed(100)
plot(exr_out)
dev.off()

#gov_approve
for_sate <- as.formula(paste0("gov_approve ~ treatment + ", paste(covariates, collapse = "+")))
lm_sate <- lm(for_sate, data = data)
sate_est <- summary(lm_sate)$coef["treatment", c(1,2)]
sate_est

exr_out <- exr(outcome = "gov_approve", 
               treatment = "treatment", 
               covariates = c("Democrat","Republican","Independents", "sex", "White", "Black",
                              "Hispanic", "Asian", "Native_american", "Middel_eastern", "Mixed", "age",
                              "education","hawk"), 
               data = data,
               boot = TRUE,
               sate_estimate = sate_est) 

plot(exr_out)

pdf(file="../Appendix/Figures/S1/C11a.pdf")
set.seed(100)
plot(exr_out)
dev.off()


#torture
for_sate <- as.formula(paste0("torture ~ treatment + ", paste(covariates, collapse = "+")))
lm_sate <- lm(for_sate, data = data)
sate_est <- summary(lm_sate)$coef["treatment", c(1,2)]
sate_est

exr_out <- exr(outcome = "torture", 
               treatment = "treatment", 
               covariates = c("Democrat","Republican","Independents","sex", "White", "Black",
                              "Hispanic", "Asian", "Native_american", "Middel_eastern", "Mixed", "age",
                              "education","hawk"), 
               data = data,
               boot = TRUE,
               sate_estimate = sate_est) 


plot(exr_out)

pdf(file="../Appendix/Figures/S1/C11c.pdf")
set.seed(100)
plot(exr_out)
dev.off()

#Create Figure B9 (Moral Licensing)
data_new <-data %>%
  mutate(.,
         torture_2=as.numeric(torture),
         condition_2=case_when(condition=="shame"~"Shaming",
                               condition=="no_shame"~"No Shaming"))%>%
  filter(!is.na(condition_2))

summary_df<-data_new%>%
  group_by(condition_2) %>%
  do(tidy(lm_robust(torture_2 ~ 1, data = .))) %>%
  mutate(torture_2 = estimate) %>% 
  drop_na()

ggplot(summary_df, aes(condition_2, torture_2)) +
  geom_point(
    data = data_new,
    position = position_jitter(width = 0.2, height = 0.1),
    alpha = 0.2
  ) +
  geom_point(size = 3, color = "firebrick3") +
  geom_line(color = "firebrick3",linetype="dashed", group=1)+
  geom_text(aes(label= round(estimate, digits = 2)), vjust=1.2, hjust=-1, size = 4, color = "firebrick3",
            family = "Times") +
  annotate("text", x = 2.4, y =1, 
           label = "paste(italic(p), \" =0.03 \")", parse = TRUE, 
           family = "Times", size = 3, color = "darkblue") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, color = "firebrick3") +
  theme_bw() +
  #scale_y_continuous(breaks = seq(0, 1, length.out = 5)) +
  coord_cartesian(ylim = c(1,7)) +
  theme(axis.title.x = element_blank()) +
  ylab("Opposition to government use of torture")

ggsave("../Appendix/Figures/S1/B9.pdf", width = 4, height = 6)

